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We formulate a fully microscopic approach to large-scale nuclear dyuamics usiug a hyperradius 
as a collective coordinate. An adiabatic potential is defined by taking account of all possible con¬ 
figurations at a fixed hyperradius, and its hyperradius dependence plays a key role in governing the 
global nuclear motion. In order to go to larger systems beyond few-body systems, we suggest basis 
functions of a microscopic multicluster model, propose a method for calculating matrix elements of 
an adiabatic Hamiltonian with use of Fourier transforms, and test its effectiveness. 


I. INTRODUCTION 

Atomic nuclei present a unique example of self-bound, finite quantum many-body systems. They not only exhibit a 
variety of excitation modes but also decay or fission into two or a few fragments. Exploring the excitation mechanism 
based on single-particle, collective and clustering degrees of freedom is an interesting subject. Intrinsically different 
shapes such as prolate-oblate may coexist or mix at close energies, leading to the so-called large-amplitude collective 
motion [l| . Spontaneous fission and sub-barrier fusion are also typical examples of the collective motion that involves 
a large-scale change of the nuclear size ii- A fully microscopic description of their dynamics is still a long-standing 
challenging problem. 

All of the above phenomena should in principle be described starting from a Hamiltonian of the system. What is 
often performed is, however, to solve an equation of motion with some constraints or to calculate energy surfaces 
assuming different shapes in order to look for a path along which the collective motion proceeds. In the case of a 
deep sub-barrier fusion an initial fragment decomposition is maintained for the whole fusion process and the relevant 
fusion potential is calculated as a function of the relative distance of the fragments. It is hard for that approach to 
take into account couplings with configurations corresponding to a different mass distribution of the fragments. Since 
the phenomena are very complicated, those approaches sound reasonable. However, neither the geometrical shape 
nor the relative distance between the fragments is a nuclear collective coordinate in a strict sense. A question thus 
arises of whether or not we can describe the large-scale dynamics by employing a true collective coordinate. 

The purpose of this paper is to make use of a hyperradius as a collective coordinate, and to step forward for a 
consistent formulation of the large-scale dynamics together with the underlying collective potential. Most of the 
needed ingredients are available in the literature. The hyperradius is a global coordinate that measures matter size, 
and it is widely used in three-body problems il. After pioneering work with the hyperspherical approach Q, its 
extension to W-body systems has been proposed for solving various problems pTl - [l3| . A common foundation of all 
the hyperspherical approaches is that a total wave function of the system is expanded in terms of a product of the 
hyperradial and hyperangular functions. There are two types of realization for describing the hyperangular functions. 
One is to use hyperspherical harmonics [Mini, and the other called an adiabatic hyperspherical approach is to 
employ channel wave functions that are defined by diagonalizing the hyperangular part of the Hamiltonian [I10. 
The former has the advantage that the hyperspherical harmonics are well-known eigenfunctions of the angular part 
of the multi-dimensional Laplacian, but its use in a real problem is fairly complicated and so far limited to few- 
particle systems. Moreover, a convergence with that expansion is rather slow. See the first paper in Ref. [ll| for 
details on the development and difficulty. The latter is widely used in atomic and molecular physics. Since an 
adiabatic hyperspherical potential defined there reflects the large-scale change of the system, we adopt the adiabatic 
hyperspherical approach in what follows. 

The equation of motion in the hyperspherical method is the same independent of the number of particles in the 
system, which is an appealing feature of the hyperspherical approach. In spite of the various efforts, only small 
systems have so far been investigated mainly because calculating the matrix element of an adiabatic Hamiltonian is 
still not trivial and solving its eigenvalue problem is hard for general A^-body systems. Correlated Gauss functions 
(CG) are employed for studying cold atom physics and electron-positron systems in Refs. [I1,[I3|, but their application 
is limited to few-body systems with the total orbital angular momentum L = 0 and 1. The use of harmonic-oscillator 
shell-model wave functions is discussed in Ref. for calculating the matrix element needed in the hyperspherical 
approach. The oscillator basis is convenient for representing such one-centered configurations that are not highly 
excited from the ground state, but it is not flexible enough to cope with a description of large-scale change such 
as, for example, the clustering and fragmentation. We instead need basis functions of cluster type to describe such 
configurations. We attempt here an extension of microscopic multicluster wave functions Ilq used to describe 
the structure of light nuclei. In the multicluster model the intrinsic fragment wave functions are described with 
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shell-model type configurations, while the relative motion among the fragments is described with the CG [13[ 13 . We 
apply a Fourier integral for evaluating the matrix element as it is applicable for any type of many-body basis function. 

The structure of the present paper is as follows. We define in Sect. [Ill the hyperspherical coordinates and separate 
the kinetic energy of the system into hyperradial and hyperangular parts. In Sect. lIIll we define the eigenvalue problem 
of the adiabatic Hamiltonian and present the equation of motion for hyperradial functions. In Sect. 113 we discuss 
qualitative features of the adiabatic potential together with a separation of active and inactive degrees of freedom. 
In Sect. we define basis functions of the multicluster model and give a method for calculating the matrix element 
integrated over the hyperangles together with examples for the overlap and kinetic energy. In Sect. I VII we show how 
to extract the evolution of intrinsic shapes of the system as a function of the hyperradius. In Sect. IVIII we touch on 
an eigenvalue problem of the full Hamiltonian with a constraint of the mean-square matter radius in comparison with 
the present approach. Conclusions are drawn in Sect. IVIlTl 


II. HYPERSPHERICAL COORDINATES 


A 


We start from defining the hyperradius for a general case consisting of K particles. The mass of the Hh particle is 
i in units of a suitable mass m. By denoting its position coordinate by Hi, we define a set of Jacobi coordinates by 


X, = 




( 1 ) 


where Ai 2 ...i = ^.nd is the reduced mass factor, /ii = Ai 2 ...iHi+i/Ai 2 ...j+i. The square of the hyperradius 

p is defined by 


P 


2 


K-l 


( 2 ) 


which is also rewritten in several ways as 

K K K 

= y^A,{R,- = y^A,Rj - Ai2...kRI^ = - V AA,{R,-Rj)\ (3) 

where i2cm is the center-of-mass (cm) coordinate of the system, ilcm = MRijA\ 2 ...k- Note that mp^ is equal 

to the trace of the moment of inertia tensor of the system. 

It is straightforward to extend the above definition to an 7V-nucleon system. Protons and neutrons are assumed to 
have an equal mass, the nucleon mass, which is taken as m. By denoting the nucleon’s position coordinate by r^, we 
define Jacobi coordinates as 


Xi = 


ri+i 


f=i 


( 4 ) 


with Pi = + Then p^ reads 

N-l N N ^ ^ 

i=l i=l i=l j>i=i 

We often use a matrix notation. For example, x = (xi) stands for an Y — 1-dimensional column vector or an 
{N — 1) X 1 matrix, and x stands for its row vector. The p^ is simply written as a scalar product, p^ = xx. It is 
clear that p^ is equally defined by any coordinates that are related to x by an orthogonal transformation. In fact p^ 
is independent of any choice of such coordinates. 

A measure of the nuclear size, p^/N is an operator for the mean-square matter radius. Symmetric with respect to 
the nucleons’ coordinates, p is a collective coordinate that has a unit of length. The other 3N — 4 coordinates are 
hyperangle coordinates denoted by H collectively. The volume element for integration reads 


dx = dx\dx2 ■.. dxjq-i = p'^ ^dpdQ, 


( 6 ) 
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where d is the dimension of the spatial coordinates excluding the cm coordinate 

d = 3{N-l). 


( 7 ) 


It is well known that the volume of a d-dimensional hypersphere with radius p = po is given by Vd = ~ 

(poV^)‘^/r(d/2 -|- 1) with the gamma function F. Since Vd is equal to p'^~^dp f dfi, the surface area of the 
hypersphere is 


/ 


dfi 


2^ 

r(d/2)- 


( 8 ) 


The volume element in the single-particle coordinates reads dr\dr 2 ■ ■ ■ drjq = N^^^dxdRcm- 
Let us introduce dimensionless coordinates by Xi = p^i. They are subject to the constraint 
An explicit form of fl may be constructed from the TV — 1 vectors but it is not needed in what follows. It should 
be noted, however, that a variety of configurations or shapes of the nucleus correspond to different functions of Q. 
The total kinetic energy T of the A^-nucleon system, with its cm kinetic energy Tcm being subtracted, is separated 
into hyperradial (Tp) and hyperangular (Tq) parts: 



2m 


N ^2 

_ T 

t—1 


f.2 a2 

i—1 '' 


(9) 


with 

T = ^ A. 

^ 2m \dp'^ p dp) 2m^ dp dp 


( 10 ) 


The hyperangular kinetic energy Tq may be expressed as 


2mp^ 


( 11 ) 


where /C(f2)^ is the square of the CTand angular momentum. An explicit form of /C(0)^ is available in a recursive way 
together with the definition of 17 [y, [I^ . 

Suppose that the iV-nucleon system develops into K fragments or clusters each of which has Ni nucleons (X)i=i 77* = 
TV). It is convenient to divide p^ of Eq. ([S]) into two groups: 


K K 

= Pin + Prel = ^ Pi + ^ Ni{Ri — Rcmf, (12) 

i=l i=l 

where Ri is the cm coordinate of the Tth fragment and p| is its squared hyperradius, 

Ni 

pl=J2{rN,,...^.,+j-Ri)^ (iVo = 0). (13) 

The first term pf^ of Eq. m gives a measure of the sum of the squared matter radii of the fragments (each pf/Nt is the 
mean-square matter radius of the Tth fragment), while the second term p^^j is exactly the same as that of Eq. ([3]) with 
Ai = TVi, giving a measure of the spatial extension of the relative motion of the fragments. It is natural to arrange the 
coordinates into cluster-internal and cluster-relative to describe the motion of the K fragments. The cluster-internal 
coordinates, denoted by (a:i, X 2 ,..., xn-k), consist of a collection of Jacobi coordinates of each fragment, and the 
cluster-relative coordinates denoted by {x]\[-k+i,xn-k+ 2 , ■ ■ ■ ,®Af-i) are Jacobi coordinates as defined by Eq. O- 
Clearly p^ is independent of the number of fragments into which the TV-nucleon system develops. 


III. EQUATION OF MOTION IN ADIABATIC HYPERSPHERICAL EXPANSION 

To solve a Schrodinger equation for the system in the hyperspherical method, a total wave function 'I' is 
usually expanded in terms of a complete set of the hyperspherical harmonics or AT-harmonics Y\(ri): d/ = 





4 


p-(rf-i)/2 [^ Here Y\{n) is an eigenfunction of Tq labeled by A. The hyperradial 

functions Xx{p) determined from a set of coupled-channels equations. This method is successfully used in nuclear 
few-body systems However, the number of hyperspherical harmonics needed to reach a converged solution be¬ 
comes very large at large p values [2l|. Moreover, the coupling matrix elements between different are of the 

same orders of magnitude as the diagonal matrix elements especially for the Coulomb interaction, which also makes 
the convergence slow. 

In low-energy phenomena, the hyperradial motion is expected to be slow compared to the hyperangular motion. 
Thus an adiabatic potential that takes account of all possible hyperangular motion at a fixed hyperradius giv es insight 
into the dynamics of the system’s evolution (2^ . We adopt the adiabatic hyperspherical expansion method BSEMl 
used extensively in atomic and molecular physics. 

We define an adiabatic Hamiltonian ifad by 


a., = r„ + r + _ T. 


Smp^ 


Smp^ 


(14) 


Here V is the total potential energy and H = T + V is the total Hamiltonian of the system. The nucleon-nucleon 
interaction of V is assumed to be an effective interaction that contains no strong short-ranged repulsion. Assuming 
that V contains no derivative of p, we solve an eigenvalue problem of the Hermitian operator i/ad, 


H,Mp,n) = u,{p)^,{p,n), (15) 

to obtain channel wave functions $;y(/3, H) and real adiabatic potentials U^{p) that are labeled by v. The quantum 
numbers of H such as spin-parity J'^ are preserved as those of i/ad, and the antisymmetry requirement on 4' applies 
on $;,(p,O) as well. Note that p appears parametrically in Eq. (ITSl) . At fixed p, all possible couplings among 
various hyperangular configurations are taken into account to obtain Uv{p) and 4)i,(p, H). The 4>;^(p, H) form a set of 
orthonormal functions at each p, 

{'^u:{p,Vt)\i^^{p,Vt))Q,=5uy, (16) 

where (.. .)a indicates that the integration is carried out over with p being fixed. Actually, the 4>jy(p, H) also contain 
spin and isospin coordinates that have to be integrated, but they are omitted for the sake of simplicity. Apparently 
Uv{p) contain the minimum ‘centrifugal potential’, h^{d— 5){d— l)/8mp^, for A > 3 even when the eigenvalue of 
/C(H)^ vanishes. 

The Schrodinger equation, = AT, is solved by expanding T in terms of 4>i/(p, fl); 


T = p-(‘^-i)/2^/,(p)4>,(p,H). 

y 


(17) 


The normalization of T is \ fi^{p)\^dp = 1 for a bound state. The hyperradial functions fv{p) are determined 

by solving a set of coupled-channels equations, 




iM-hT. 


‘2-Puv' (p) -1" Qvv' (p) 

dp 


fw (p) = 0, 


with non-adiabatic coupling terms 


o 2 

Puy'ip) = (Ti.(p,fl)|^Ti.'(p,n))n, Quu'ip) = {<^^ip,n)\-^<P^>{p,n))Q. 


(18) 


(19) 


Equations (I15L (TTSl) . and m give a microscopic description of the large-scale dynamics. 

A unique advantage of the adiabatic hyp ersp herical approach is that both lower and upper bounds to the exact 
lowest energy of H are readily obtained [231124| . As shown in Appendix I3 we have 


Tbi/(p) = 0 . 


( 20 ) 


Its differentiation with respect to p leads to 

Qw^P) Y Ql'V^P') 0 : 


( 21 ) 


where 


Qvu{p) = {^^^u{p,Yt)\—^y{p,n))n 


( 22 ) 
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is non-negative, and consequently Qvv{p) < 0. The potential defined by 

Wi/(p) = Uv{p) — 

always satisfies W^{p) > Uv{p). The lowest eigenvalue E obtained by truncating Eq. (fTSl) to a single-channel equation 
with the lowest adiabatic potential Uo{p) or Wo(p) gives a lower or upper bound to the exact lowest energy of H. See 
Appendix for details. Convergence of the solution of Eq. (fTSl) is checked by increasing the number u of channels. 

A time-dependent Schrddinger equation is convenient for studying how final configurations in e.g. few-body decay 
and sub-barrier fusion evolve from their initial states. The wave function at time t is assumed as 

^ 12 ). ( 24 ) 

1/ 


Once fv{p,0) are given, fu{p,i) for t > 0 are determined from the equation 


d 


2m dp'^ 


uAp) 






2Tbi/'(p)^—I" Quv'{p) 

dp 


U'iP,t)- 


(25) 


IV. HYPERRADIUS DEPENDENCE OF ADIABATIC POTENTIAL 

The p-dependence of Uu{p) or MA(p) governs how the nucleus responds to its change of size. The kinetic energy and 
the Coulomb potential respectively give 1/p^ and 1 /p contributions to Ui,{p) at large p values. Short-range pairwise 
nuclear interactions give a p~'^{n > 3) contribution |25l| . Let us focus on the lowest adiabatic potential with the same 
spin-parity J'^ as that of the ground state. Uo(p) has a minimum at p ~ pmin corresponding to the matter size of 
the ground state. As p decreases from pmin, Uoip) rises because of a loss of nuclear potential energy as well as an 
increase in the kinetic energy. As p increases from pminj various configurations contribute to determining $o(p, 12). 
Here, deformations, shell effects, couplings with different modes and so on participate in determining the adiabatic 
potential. Uo{p) reaches a peak at some p value or may even have a couple of local peaks at different p values. As p 
increases further, t/o(p) approaches the lowest decay threshold of the nucleus. 

The above global feature of the adiabatic potential well corresponds to the decomposition of p^ in conformity with 
a formation of fragments or clusters. As shown in Eq. (1121) . the different decomposition of the fragments can be 
treated on an equal footing in the hyperspherical approach, which makes it possible to assess what configurations 
play an important role in determining the adiabatic potential. If one instead calculates a sort of adiabatic potential 
or potential energy surface as a function of the relative distance between two fragments, there is no way to compare 
such potentials for different fragment decompositions because their relative distances have a different meaning. 

What fragment decompositions or configurations contribute to the lowest adiabatic potential clearly depends on p. 
The expectation value of iL is a major contribution to the adiabatic potential (see Eq. (03 1. We rewrite H according 
to the fragment decomposition: 


K 

H = Y,{T, + V,)+T,,i + V,eu (26) 

i=l 

where Ti + Vi is the intrinsic Hamiltonian of the Ah fragment, Tj-ei the kinetic energy of the relative motion among the 
fragments, and I4ei denotes the potential energies acting between the nucleons belonging to the different fragments. 
Vrei depends on both cluster-internal and cluster-relative coordinates, thus causing a coupling of the relative motion 
among the fragments with their intrinsic motion. When p^ei is so large compared to pin that the nucleon-nucleon 
interactions of Kei can be neglected and only the leading term of the Coulomb potentials of Kei is retained, l^ei 
reduces to 


K 

Hrel —t 

j>i=l 


I ~ ^3 I 


— C(Hrel), 
Pre\ 


(27) 


where ZiC is the charge of the ith fragment and prei and flrei are the hyperradius and hyperangles constructed from the 
cluster-relative coordinates {x^-k+i, xn-k+ 2 , ■ ■ ■, xn-i)- With increasing p the intrinsic motion of each fragment 
is stabilized toward its own ground state, while the configurations responsible for the relative motion are decoupled 
from the intrinsic motion. Both the coupling and decoupling of various degrees of freedom are naturally taken into 
account in the hyperspherical approach. 
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FIG. 1: The 10 lowest adiabatic potential curves of the three-a system with J'^ = O'*" that are taken from Ref. [^. The 
hyperradius denoted R here is defined by = %/3 — Rcmf', and the energy is measured from the three-a threshold. 

The solid line denotes the adiabatic potential dominated by the ®Be+a channel, while the dashed lines denote the potentials 
dominated by the three-a continuum channel. Panel (b) is an enlarged view of the potentials at large R where a number of 
sharp avoided crossings successively appear. 


When there are several thresholds corresponding to different fragment decompositions, avoided crossings of the 
adiabatic potential energy curves may occur. As an example, we show the case of that is described with a cluster 
model of three a-particles [2^. The eigenvalue problem ((TS|) for i/ad is solved accurately, and an analysis of the 
adiabatic potentials clarifies how the contributions of the hyperangular kinetic energy, the nuclear potential and the 
Coulomb potential change as a function of p. Figure[T] taken from Fig. 2 of Ref. [l^, displays the 10 lowest adiabatic 
potential curves for = 0"*". The lowest potential Uq{R) has a minimum at i? « 3.5fm, which is deep enough to 
support a bound state, that is, the ground state of Furthermore, the lowest potential reaches a broad peak around 
12 fm, corresponding to the second 0+ state of the Hoyle resonance state. The adiabatic potential indicated by 
the solid line is dominated by the two-body ®Be-|-a state and approaches the ®Be-|-a threshold at large R, while the 
other potentials indicated by the dashed lines are all dominated by the three-a continuum states. As seen in Fig. [T] 
(b), an avoided crossing begins to occur at i? ~ 140 fm, which is because the three-a continuum state comes down 
closely to the two-body ®Be-|-a state. Since the avoided crossing actually occurs within a small range of R, it may 
be hard to see it in the figure. Refer to Fig. 3 of Ref. [2^ to confirm the crossing clearly. Since the ®Be-|-a threshold 
is higher than the three-a threshold, a number of avoided crossings successively appear below the adiabatic potential 
indicated by the solid line. As is well known, the non-adiabatic coupling terms (11911 may be singular especially when 
the avoided crossing is sharp, narnely it occurs within a small range of p. In that case, a diabatic procedure is proposed 
for accurately solving Eq. (fTSl) [a . The slow variable discretization method combined with a complex absorbing 

potential makes it possible to solve Eq. (fT^ and to reproduce the energy and width of the Hoyle resonance in good 
agreement with experiment [^ . 

Let us speculate concerning the adiabatic potential curves of ^g|Cf that are crucially important for determining its 
decay mode. The ground state of ^®^Cf decays mostly by an a-particle emission. The rest is a spontaneous fission 
(SF), emitting 3.7 neutrons on average. To make things simple, we approximate the SF as occurring through a single 
channel of ^ 54 X 0 -I- ^^iRu -I- 4n. The two decay modes contain different numbers of fragments, two in a -I- ^|gCm 
and six in the SF, but the hyperspherical approach can treat both in a unified way. The threshold of a-|-^"^®Cm is 
6.2 MeV below the ground state of ^®^Cf, whereas that of the SF is 200.4 MeV lower than the ground state. See the 
schematic diagram of Fig. [2J The lowest adiabatic potential Uq{p) approaches the SF threshold at large p. Above 
that threshold many U^{p) curves, not drawn in Fig. [2J show up corresponding to the continuum states of the SF 
mode. A unique U^{p) with the two-body a-|-^'^®Cm character appears high above the SF threshold. When moving 
inward from this asymptotic region, the Coulomb potential ( 1271 ) produces a distinct difference between the two decay 
modes. The charge factor Z\Z 2 of the SF mode is more than ten times larger than that of the a channel. Thus those 
Uv{p) curves that are dominantly contributed by the SF configurations rise up rapidly, while the U,^{p) curve of the 
a channel increases much more slowly. At the avoided crossing point pacj the lowest curve Uo{p) comes very close to 
that of the a curve, and for p < pac the a channel makes a dominant contribution to Uo{p). With further decrease of 
p many different configurations begin to mix due to an increasing role of the nuclear interaction V. The Uo{p) reaches 
a barrier top around some point and reaches its minimum at pmin corresponding to the matter radius of the ground 
state of ^®^Cf. Though much more complicated than the case, the gross feature of the adiabatic potential curves 
of ^||Cf should have some similarity to those of and the decay branch of ^||Cf will be determined by solving 
Eq. ( 1151 ) . 
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FIG. 2: A schematic diagram of the adiabatic potential curves of ^HCf as a function of hyperradius p. Energy is measured 
from the ground state of ^ggCf. The potential dominated by the a + ^gfCm channel (solid line) goes to —6.2 MeV at large p 
values, while that by the ^ 54 Xe + ^ffRu + dn channel (dashed line) to —200.4 MeV. An avoided crossing occurs at p « pac. The 
lowest adiabatic potential Uo{p) changes its dominant character from a + ^gfCm to the ^I^XeT ^ 4 |Ru + 4n channel around pac. 


V. MULTICLUSTER APPROXIMATION AND INTEGRATION OVER HYPERANGLES 

Solving Eq. (113 is of vital importance in the adiabatic hyperspherical approach. Its accurate solution is obviously 
very hard except for few-body system. The difficulty is enhanced by the fact that the matrix element has to be 
calculated by integrating over U only. Some efforts have been made for extending to larger systems EElil- We 
take up this problem assuming the use of many-body wave functions that contain all the coordinates. 

Before discussing the eigenvalue problem (1151) , we note that a usual approach defines an adiabatic potential barrier 
or energy surface at a given ‘collective’ coordinate by searching for a minimum of V for various parameters that 
characterize the nuclear density or shape [2^. This makes sense in that V is a major part of H^d, and because, 
since U is a function of p and U, its minimum gives information on the most important U values contributing to the 
lowest adiabatic potential. As mentioned before, the adiabatic hyperspherical approach can go beyond that by taking 
account of various couplings with different degrees of freedom. 

Let us assume that the channel wave function d)i,(p, fl) at a given p is expanded in terms of suitable basis functions 

$j.(p, U) = ^ C^i{p)(l)i{x). (28) 

i 

Equation ITSl) is then reduced to the following generalized eigenvalue equation for determining the coefficients C,^i (p) 
and the adiabatic potential U,^{p): 

Y^imjip) - U,{p)B.M]C.j{p) = 0, (29) 

j 

where T-Lij (p) and Bij (p) are adiabatic Hamiltonian and overlap matrices defined by 

'Bijip) = {(j)i{x)\Hi,d\(t)]{x))n, Bij{p) = {(i)i{x)\(j)j{x))a- (30) 

We include only those basis functions that give a c-number p^ for the expectation value of the squared hyperradius 
operator xx'. 

{(i3^{x)\xx\(i3,{x)) 2 /oi') 

We face two problems. One is what basis functions we use for (j>i{x). The other is how to calculate the matrix 
element in Eq. (1501) . The first one is crucially important for assessing the quality of $j,(p, 0) and Though it is 

difficult to give a general answer, our ansatz is to employ a microscopic multicluster approximation pisl . This is 
because, as mentioned in Sects. HI and IIVI the structure change we are interested in includes a variety of configurations 
ranging from one-centered shell-model wave functions to those with a few fragments or subsystems. A general form 
of the multicluster wave function containing K fragments reads 

(^(^■^(a;) = A{^i{zi)^2iZ2) ■ ■ • 'i'KizK)xixN -K+lj • ■ • ; ®Af-l)}, 


(32) 
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where A is an antisymmetrizer, 'Pi an antisymmetrized intrinsic state of the ith fragment containing JVi nucleons 
and X is the relative motion function for the fragments. The cluster-internal coordinates {xi,X 2 , ■ ■ ■ are 

abbreviated as (zi, 22 ,..., zk), where e.g. Zi stands for the first — 1 Jacobi coordinates {xi,X 2 , ■ ■ ■, ajATj-i). The 
spin-isospin coordinates are again suppressed. In general may represent not only the ground state of the fragment 
but also its excited state. The quantum numbers for characterizing are omitted. The coupling of the angular 
momenta of the d'^s and y to a total angular momentum JM is implicitly understood in Eq. ([32]). We presume (j)i{x) 
to belong to the space spanned by 

+ {(j)^'^\x)} + ... (33) 

Note that any states in {(j)^^^{x)} are in general nonorthogonal to each other even when they belong to different 
K subspaces. The questions of what intrinsic states of the fragments are important and what K subspaces have 
to be included depend on a given system and energy range of interest. To proceed further, we assume that is 
approximated by harmonic-oscillator shell-model wave functions, while x is described well with a superposition of 
Gauss functions [13, [H, as developed in few-body problems. 

We have to calculate a matrix element for some operator 0{x), 

O{po) = {(j)i{x)\0{x)\(j)j{x))a^p^p^, (34) 

by integrating over PI at fixed p, say po- Tii6 calculation of the matrix element of Tp in iJad can be aided with use of 
the identity 

-*,W = -(^x..—)i,,W. (35) 

^ ^ i—1 

See Appendix iBl for an example. In calculating the matrix element of H, the cluster-intrinsic term + ^i) (see 

Eq. (|26ll l may be replaced by 

H(j)^^\x) f {x) 

-b A{{T,el + Vre\)'Plizi)'p2{Z2) ' ' ' 'P K {ZK)xix N-K+1, ■ ■ ■ , (36) 

using the observed energy Ei of 4'^. This approximation looks reasonable and practically useful because any nuclear 
interaction can not satisfactorily reproduce the saturation property of nuclear binding energies despite the fact that 
reproducing the threshold energy for the fragment decomposition is important in the present approach. 

The second problem has so far been examined using integral transform techniques [l^ . 1^ . We use a <5 function 
technique as in Ref. [l^. Using the expression for Dirac <5 function 

S{p-po) = - / AAi-pVpDdco, (37) 

7rpo J-oo 

we can express O(po) as a Fourier transform of ^Po (oj): 

O{po) = - [ A^Fp,{io)duj, (38) 

Fp„{oj)= y'e-“«((/..(pol))*0(poO</'.(Pol)dS. (39) 

Note that Xi is changed to po^i with a dimensionless variable ^i. In Eq. (15^1) stands for d^id ^2 ■ ■ ■ d^N-i, where 
the integration range of each ^i covers the whole three-dimensional space. Since e' e results in a 

simple modification of the basis function, Fpg (uj) can be calculated with a technique developed in microscopic cluster 
models [HIsS]. 

In some cases the Fourier integral (1551) can easily be obtained by Cauchy’s integral formula that reduces to a 
residue calculation. Whether or not we have a practical means for evaluating Eq. (I34|) for a general case depends on 
how fast and accurately the Fourier integral is computed. For this aim we test the Whittaker cardinal series or the 
Whittaker-Shannon interpolation formula [s^: 

CO 

Ppo{^)= Fpoi‘^n)sincf{uj-Wr,), 


n— — oo 


(40) 
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where sinca; is the sine function, sincc/a:, and a;„ = nh (n = 0, ±1, ±2,...) is the grid of the sampling points. The 
series (HOI) is known to converge if Fpg{u!) is a band-limited function. Because sinenTr = Sn^, the series is exact at 
all the sampling points. It is in fact an expansion in terms of orthogonal functions {sinc^(w — w„)} that have the 
properties: 


/ OO 

sinc^(a; — u!n)duj = h, 

-CO 

/ CO 

sincf (w - ujjn) sincf (w - = h6m,n, 

-CO 


e*‘^sinc^(a; — ujn)duj = (^ < 1). 


(41) 


The third equation called the Dirichlet integral leads to an approximation for O{po): 


h _ 

TT ^ 
n— — M 


(42) 


which is nothing but a trapezoidal rule for the integration. This result is due to the fact that the Fourier transform 
of the sine function is the rectangular function and vice versa. To determine M, we need to know how fast Fp^ (w) 
decreases as a function of w. The mesh size h{h < n) is determined by examining how accurate the expansion is at, 
e.g. uj = {n + ^)h, the midpoint of ujn and tOn+i- 

Other interpolations, e.g. a cubic spline interpolation may also be worthwhile testing because it leads to a simple 
expression for Eq. (13811 and in addition the mesh size can be taken as piecewise variable. Once dFpg{(jj)/duj values 
at both boundaries of the interpolation are calculated, we can completely fix the interpolating function of the cubic 
spline. 

Since w-dependence of Fpg{u}) is of practical importance, we examine it for the diagonal matrix elements = 

4>j{x)) of 0{x) = 1 and Tn in a very schematic model. As the model, we employ CO ignoring the antisymmetry 
requirement of the wave function and focus only on its spatial part. See Appendix [B] for some basic matrix elements 
with the CG. For a spherical CG, exp(—^alAa;), the positive-definite symmetric matrix A is set to TrA“^ = 
because of Eqs. dSO and (IB2I) . We may choose A to be diagonal, A = {aiSij), as far as the diagonal matrix element 
of 0{pq) is concerned. 

Our first choice for A is a uniform nuclear expansion, at = a, leading to a hyperscalar Gaussian, exp(—= 
exp(—iap^). This function is totally symmetric and fl-independent. By taking a as {N — l)/a = |pq, the overlap 
matrix element is (see Eq. (IB3I1 1 


/ /p 'lAT-l \ 3/2 / \ d/2 

Glearly |F’pg(u;)| becomes very small if ui is significantly larger than d/2. The Fourier transform (155)) can be rigorously 
computed in this case. If d/2 is an integer, Fpg{uj) has a pole of order d/2 at w = id/2, so that the integral is reduced 
to a residue calculation, yielding 


O{po) 


-d/2 

r{d/2) 


(44) 


Even when d/2 is a half integer, we can derive the above result as follows. By the change of the integration variable, 
d/2 + ioj = —t, O{po) is reduced to 


O{po) 



p — d/2—ico 
' —d/2-\-ico 




(45) 


By changing the integration path to the Hankel contour and using Hankel’s integral representation and Euler’s 
reflection formula for the gamma function, we find the above integral to be 27r/iF((i/2). The result (H5)) is in fact 
trivial thanks to Eq. ([8|) if we note that the hyperscalar Gaussian at p = po is exp(—^apg) = and hence O(po) 

must be f dfi. We note that O(po) for 0(x) = Tn vanishes because the hyperscalar Gaussian is 11-independent 
and, when acted on by Tq, vanishes. This is also confirmed by using Eq. (IBG) . 
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FIG. 3: The ratio r(uj) = |^po(^)/-^po( 0)I for ths overlap corresponding to N nucleons’ symmetric fission as a function of uj. 
po is set to rj'/NR^uisiN): Solid, dashed, and dotted lines correspond to rj — 1,3, and 5, respectively. See text for details. 



FIG. 4: The same as Fig. [3] but for the hyperangular kinetic energy. 


The next example is a ‘symmetric fission’, that is, the nucleus fissions into two identical fragments with mass 
number N/2 and only the relative motion between them expands with increasing po- Let i?rms(./V) denote the root- 
mean-square radius of a nucleus with mass number N, and set it equal to (tq = 1.1 fm). The matrix A 

for the symmetric fission is chosen as ai = 02 = ... = aN-2 = cl, and a and a^-i are determined by the condition 


N-2 1 2 , N/2 - 1 2 , 

-+-= oPI — -= 7:p} 

a aN-i 'J a 6 


(46) 


where pf is fixed to N/2R,-„^s{N/2). The mass number N is changed to 4, 40, and 240, and for each N pQ is 
taken as po = rj'/NR^^siN) {rj = 1,3,5). Figure [3] displays log^g ?'(w) for the overlap, where r(a;) = |F’pQ(a;)/FpQ(0)|. 
In the case of = 4, the fall-off of r(w) is slow with increasing w and 77 . For example, logj^g r{uj) at a; = 5000 is 
— 13.7, —11.4, —10.1 for 77 = 1,3,5, respectively. For N = 40, r{oj) rapidly drops to 10“^® as a function of w, but 
its decrease becomes slower for 77 = 5. This behavior is also valid for N = 240, and the decrease in uj becomes even 
slower with increasing rj. As shown in Fig. 21 the ratio logj^g r(uj) for Tq is very similar to that of the overlap. 


VI. EVOLUTION OF INTRINSIC SHAPES 

It is interesting to know how an intrinsic shape of the nucleus changes with increasing p. When a decay or an SF 
is considered as a tunneling through a barrier, the shape will give insight into where the fragments are formed and 
how they evolve during the passing through the barrier. The barrier is conventionally calculated by assuming some 
density distribution constrained with shape or deformation parameters such as quadrupole and octupole [3, [1^ . Such 
deformations are not observable, however. Our view is to reverse this approach. Since the nucleus should in principle 
preserve the total angular momentum, it is not trivial to imagine the intrinsic shape in the space-fixed frame. For 
example, any state with L = 0 is spherical in that frame, but it can happen that such state is intrinsically deformed 
and rotates. As shown in Ref. [s^, the intrinsic two-a structure of the rotational state of ®Be emerges from the wave 
function obtained by a quantum Monte Carlo calculation. 

Following the procedure of Ref. [s^, we can get the intrinsic density or deformation indicated by e.g. the lowest 
channel wave function <i)o(p, Id). Since 4)o(p, fl) is normalized as in Eq. (flBl) . its square, Fp(fl) = |4>o(p, f^)P, gives the 
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probability density as a function of at a given p. First, we generate many sampling points fli, 172 , • ■ ■, according 
to the distribution of -Pp(17) using the Metropolis-Hastings algorithm. Secondly, we define a body-fixed intrinsic 
frame for each ^lj = (£i, ■ • ■ > follows. By using Eq. ([4]) together with ilcm = 0, Jacobi coordinates 

xl = p^l (i = 1, 2,..., iV — 1) specify the positions of N nucleons (r^, ..., in the space-fixed frame. From 

these position vectors, we calculate the moment of inertia tensor 


— 


N 

2=1 


(47) 


where {a = x,y,z) is the Cartesian component of r^. Diagonalizing the 3x3 symmetric matrix determines 
the principal moments of inertia, which define the axes of the intrinsic frame. For example, the axis is called x', y', z' 
in increasing order of the principal moment of inertia. The direction of the axis also has to be chosen consistently. 
By reading as (r'j,r' 2 , ■ ■. , 1 ''^^) in reference to the intrinsic frame, we obtain the desired position 

coordinates of N nucleons in the intrinsic frame. Finally, accumulating these position coordinates over j = 1, 2,..., M 
leads to the intrinsic single-particle density at p. Once the intrinsic density is obtained, it is easy to extract multipole 
deformations. 


VII. EIGENVALUE PROBLEM OF HAMILTONIAN WITH RADIUS CONSTRAINT 

It looks as though the adiabatic hyperspherical approach has some relationship to an eigenvalue problem of the 
Hamiltonian with a constraint . Let us attempt to find a solution of the Schrodinger equation by first constraining 
the expectation value of the squared hyperradius to a fixed p^ value. Suppose that the solution is expanded in terms 
of some basis functions: 

^ G^,i{p)(j)^{x). (48) 

i 

The constraint (I5T]) is not necessarily imposed on (j)i{x) itself, but we demand the solution we are looking for to satisfy 
the condition 

p^{x)\xx\'ilp^{x)) ^ 2 /.O') 

Here k is a label related to the Lagrange multiplier. The unknown coefficients G^iip) and the energy eigenvalue E^^p) 
are determined from the following equation 

[Hij - K{Qij - p'^Bij) - E^{p)Bij\ G^j{p) = 0, (50) 

j 

where H ^ and Q are matrices defined by 


Hij = {(j)i{x)\H\(j)j{x)), Bij = {(j)i{x)\(j)j{x)), Qij = {(l)i{x)\xx\(j)j{x)). (51) 

Unlike Eq. (1501) . the above matrices are obtained by integrating over the whole coordinates. To determine the 
coefficients G^j{p) from Eq. (ISOl) . the value of k has to be given. Actually k should be such that both Eqs. (l4^ and 
(ism) are simultaneously met. Apparently 4tp,^,(a;) and pK,{x) are not orthogonal to each other even for p' = p. 

The next step is to use the generator coordinate method in which a solution 'k for the Schrodinger equation is 
assumed as 

= y] /c„(p)4'p^(a;)dp. (52) 

The coefficients Gk{p) are determined from the Hill-Wheeler equation 

y] f{'l'p>^ix)\H-E\^p,ix))G,ip)dp = 0, (53) 

K, 

which should be satisfied for any p' and k! values. An approximate solution to the Hill-Wheeler equation gives an 
upper bound to the ground-state energy. Note that the adiabatic hyperspherical approach gives both lower and upper 
bounds as discussed in Sect. uni 
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We refer to two interesting calculations with a constraint in comparison to the adiabatic hyperspherical approach. 
One is a Hartree-Fock-Bogoliubov calculation performed by constraining the mean-square radius, ^ 

study how self-conjugate nuclei fragment into a clusters [^. As Eq. ([5]) indicates, this constraint is equivalent to 
that of provided the contribution of to the squared radius remains a constant. The treatment of the cm 
motion in Ref. does not satisfy this condition as usual in a mean-filed model. It would be a challenge for the 
mean-held approximation to cope with such diverse structure at large distances that is composed of different numbers 
of fragments. What should be further pursued at this moment is to establish the essential relationship between the 
adiabatic hyperspherical approach and ‘beyond mean-held’ calculations or conhguration interaction calculations that 
constrain the mean-square matter radius. 

Another is a simultaneous study of both a-|-®He reactions and the structure change of ^°Be in a microscopic 
a-|-a-|-n-|-n model , in which a distance parameter between the two a-clusters is constrained. Since the motion of the 
two neutrons is restricted to either molecular or atomic orbits around the a-clusters, the main conhgurations included 
are a-|-®He and ®He-|-^He two-body types. The adiabatic energy surfaces are calculated within that approximation. 
An avoided crossing is treated by the generator coordinate method. As noted in Sect. lIVi the relative distance of the 
fragments is not a collective coordinate. If one constrains as the generator coordinate, it would be possible in the 
same four-body model to take account of possible couplings with the ®Be-|-n channel that is the lowest threshold of 
^°Be as well as the three- and four-body channels, ®Be-|-n -I- n and a + a + n + n, that are open in the energy region 
treated in Ref. 


VIII. CONCLUSION 

Stressing that the hyperradius is a collective coordinate, we have formulated a fully microscopic adiabatic hyper¬ 
spherical approach to large-scale nuclear dynamics. The equation of motion for hyperradial functions is universal, 
independent of the number of nucleons, and enables one to consistently treat the dynamics from confined nuclear 
motion to relative motion among fragments in their asymptotic region. It is possible to describe in a unihed way 
cases where the nucleus fragments into several channels. No spurious center-of-mass motion appears and couplings 
with different degrees of freedom can naturally be taken into account. These properties are due to the fact that 
both the squared hyperradius and the kinetic energy are flexibly decomposed into cluster-internal and cluster-relative 
quantities responding to the fragment formation. 

The adiabatic potential as a function of the hyperradius plays a key role in the present approach. It is unambigu¬ 
ously defined solely by the Hamiltonian of the system, and there is no need to assume specific geometrical shapes 
or deformations to compute it. Conversely the shape or intrinsic density, if necessary, comes out after the adiabatic 
potential is obtained or the equation of motion for the hyperradial functions is solved. The calculation of the adia¬ 
batic potential involves the integration over all the coordinates but the hyperradius. Expecting that a microscopic 
multicluster model is a promising candidate for applying the present approach to larger systems, we have discussed 
the use of Eourier transforms for evaluating the matrix elements needed to obtain the adiabatic potential. The ma¬ 
trix elements can be obtained in exactly the same way as the usual matrix elements needed in nuclear many-body 
calculations. A merit of the Fourier transform technique is its simplicity, and test calculations indicate that accurate 
evaluations of the matrix elements are feasible. 

Although the calculation of the adiabatic potential still requires much computer time for large systems, a real 
challenge is whether we can provide large enough basis functions to cover important configurations for fixed p. Further 
developments are certainly indispensable for a microscopic, realistic description of large-scale nuclear dynamics. 
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Appendix A: Lower and upper bounds 

In this appendix we rigorously prove Eq. (I20|) and show that both lower and upper bounds to the ground-state 
energy are respectively obtained by solving single-channel equations. 
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The ground-state wave function may be expressed in the hyperspherical coordinates as 

= (Al) 


with the normalization condition 


pOO 

/ I/(P)P 

^0 


dp=l, ($(p,fl)|$(p,n))n = 1 . 


The hyperradial function f{p) has to vanish at p = 0. The ground-state energy reads 


& 




d'^f{p) , ^df{p) 


dp'^ 


+ 2^p(p) + /(p)g(p)jdp + \fip)\^u{p)dp, 


where 


r\ 

P{p) = mp,n)\—^p,Lu))n, Qip) = (<i>(p,fl)|^d>(p,f2))n, 


U{p) = {^pM)\H,<impM))n■ 

F^:o'^l the normalization condition of $(p, fl), we obtain 


(A2) 


(A3) 


(A4) 


^mp, fl)|$(p, fl))n = P{pr + Pip) = 0. (A5) 

Thus Pip) must be pure imaginary or zero. If Pip) is not zero but pure imaginary, f*{p)df{p)/dp in Eq. (IA3I) must 
also be pure imaginary because Aexact is real. With f{p) = g{p) + ihip), where gip) and h{p) are real functions, 
f*ip)dfip)/dp reads 


^ 1 A{g(p )2 i{g(^p)h'{p) - h{p)g'{p)}, (A6) 

which leads to dlgip)"^ + /i(p)^}/dp = 0. Thus gip)^ -I- h{p)‘^ is a constant, and it must be zero because of /(O) = 0. 
Namely, f{p) vanishes identically, which can not be accepted. Using P{p) = 0 in Eq. (jA3l) leads to 

= f f-M {-~ + ir(rt) imp (AT, 


with 


Wip) = Uip)-^Qip). (A8) 

Suppose that for $(p, U) we take the d>o(p, ^) that gives the lowest adiabatic potential. The corresponding quantities 
Qip) and U{p) in Eq. (IA4I) are denoted by Qoip) and Uoip), respectively. It follows from the Ritz variational principle 
that 


inexact < I rip) + Woip)j fip)dp. (A9) 

If f{p) is chosen to be the solution of the equation (the adiabatic approximation), 

with the lowest eigenvalue Ejj, E\j turns out to be an upper bound of Eexact^ Inexact < E\j. Differentiating Pip) = 0 
with respect to p leads to 


d d 

+ Qip) = 0- 


(All) 
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Equation (ITtI) for ^exact is rccast to 


-E'fixact — 


np) 


cP 


+ Uip) 


2m dp^ ' “ \/-y ' 2m^ dp"dp 


f{p)d-P- 


Since the last term in the square brackets is non-negative, we obtain 

+c'(o) I nodf^ 


inexact r{p)[- 


(A12) 


(A13) 


By using the inequality U{p) > Uo{p) and choosing f{p) to be the solution of the equation (the Born-Oppenheimer 
approximation), 


(-|^^ + t/o(p)) fip) = E^fip), (A14) 

with the lowest eigenvalue Al, we obtain a lower bound of Aexact as Aexact > Al- 

If we calculate the expectation value of H for the wave function 4^ = f{p)^^{p,n) with the z/th channel 

wave function, we confirm Eq. (EHl) using the same argument as above. 


Appendix B: Matrix elements with correlated Gaussians 

In this appendix we calculate Fp^{uj), Eq. dSi), using as 4>i{x) the generating function g{s]A,x) [l3, S) HI) 113 of 
the CG: 


g(s; A, a;) = exp(—iiAa: -I- sa:), (Bl) 

where A is an (N — 1) x {N — 1) symmetric, positive-definite matrix and s = (si) is an (N — l)-dimensional column 
vector to describe motion with non-zero orbital angular momentum. They are both parameters that characterize the 
CG. The constraint (1^ reads 

^TrA-^ -h sA-^s Ri pI (B2) 

Note that for the special case that A is diagonal, A = (aiSij), g{s\A,x) reduces to a product of Gaussian wave 
packets: g{s;A,x) = n„--T^ exp[— ^a,:(a^?; — + ^Uisf], We present formulas for (a;) calculated between 

g{s-,A,x) and g{s';A\x). See Ref. [l^ for details. The case with s = s' = 0 is given in Ref. pOl| . 


Overlap 

The function Fpg{u!) for 0{x) = 1 is given by 


with 


Fpo(^) 


(2^)^-i 

deti? 


3/2 


— ■kvB 


B = pI{A +A')+ 2iujl, v = pq{s + s'). 


(B3) 


(B4) 


Here I is the {N — 1) x (iV — 1) identity matrix. Since A -|- A' can be diagonalized by an orthogonal matrix, the 
matrix B can be diagonalized as well. 


Kinetic energy 

To calculate the matrix element for 0{x) = Tq = T — Tp, we use the following relation [i^ : 

d I 

— 5 (s; A, x) = -{-xAx -H sx)g{s] A, a;), 

9^ I r 1 

7 ^g(s; A, a;) = — (xAx)'^ — xAx + (sa;)^ — 2{xAx)sx g{s-, A, a;), 
dp^ L J 

fP r 1 

Tg(s-,A,x) = - — ATvA + ss — 2sAx + xA^x g(s]A^x). 

2m 1 J 


(B5) 









15 


Combining these results, we obtain 


Tag{s-,A,x) = 


2mp^ 


Sp^TrA — p^as + 2p^sAx + {d — l)sa; 


+ {sx)'^ — p^xA^x — dxAx — 2(xAx)sx + (xAx)'^ g{s-, A, x). 
Here d is defined in Eq. ©• The change of variables, x —>• yields Fp^iu) as 






3poTrA - pIss + 2plsA$ + {d - l)pos$, 


+ - dpliA^ - 2pi{iAi)s^ + ptimf] d^. 

This integral can be performed analytically. In the case of s = s' = 0, we obtain 


P 1 1- 

2mpl 


(27r) 


7V-1 \ 3/2 


detH 


STrp^A - 3TrH-^(p^A)^ - 3dTrH-Vo^ 


+ 15(TrH-Vo^)^ - 12M2{B,pIA) 

Here use is made of the formula 


e-^^^^{iPxYdx = 


(27r) 


N-l \ 3/2 


deti? 


15(TrH-^P)"’ - 12M2{B,P) 


M2{B,P) = 


detB 


where 

M2{B,P) 

is defined by 



Bii 

Bi2 

. Hi 

N-l 

N-l 

U.i 

P^2 ■ 

. P^ 

N-l 

E 





j>i=l 

Pp 

Pp . 

■ Pi 

N-l 


U 7 V-II Bj\f-12 ■ 

■ Bm- 

-17V-1 


(B6) 


(B7) 


(B 8 ) 


(B9) 


(BIO) 


Potential energy 

The matrix element for 0{x) = V is conveniently calculated by expressing the distance vector of two nucleons as a 
combination of Jacobi coordinates 


(Bll) 

where C, is an {N — l)-dimensional column vector determined by i and j. For a Gauss potential, 0{x) = = 

^-PxCCx^ E'pj,(w) reduces to that of the overlap. For s = s' = 0, we obtain 


Fpoi^) = 


(27r) 


N-l 


3/2 


(27r) 


N-l 


3/2 


det(H 


2 r2poCC). 


(1 + 2t2poC 5-^C) detH 


(B12) 


In the last step Sherman-Morrison formula, det(i3 + c^C) = (1 + cCB~^Q detB, is used, where c is a constant. 

By comparing this result with Eq. (IB3I) and by noting that B~^ is constrained by the condition (IB2I) . we expect 
that the contribution of the potential energy V to the adiabatic potential behaves as {p/t)~^ for large p values. The 
p~^ dependence was found for a three-body system [^ . but our result suggests that it is valid for many-body systems 
as well. 

As an important application of Eq. (IB12I) . we calculate the matrix element for the Coulomb potential, 0{x) = 
1/ki - ’’jl- Using 


= — / 


'TT Jo 


(B13) 
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and the integral 


'OO 

(1 + ar^y^/^dr 


1 

7 ^’ 


(B14) 


0 


we obtain the matrix element for the Coulomb potential as 



(B15) 


As expected, the inverse p-dependence appears naturally. 
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